***************************************************************************
***************************************************************************
**Replication files for Enns, Moehlecke and Wlezien: "Detecting True 
**Relationships in Time Series Data with Different Orders of Integration"**
**Table A-2, appendix******************************************************
**Last run: Feb 6th, 2021************************************************** 
***************************************************************************
***************************************************************************



****************************************
**Simulation Code to Replicate Table A2
****************************************

**********************************
*Table A-2 layout
**********************************
*generate excel file to store Table 1 output
putexcel set tablea2.xlsx, sheet(sheet1) replace
*add labels to table

putexcel A4 = "$\hat{\alpha}_1$"
putexcel A5 = "$\hat{\beta}_1$"
putexcel A6 = "$\hat{\beta}_2$"


putexcel B1 = "T = 50"
putexcel B2 = "$\rho = 0.2$"
putexcel B3 = "coef"
putexcel C3 = "%"
putexcel D2 = "$\rho = 0.5$"
putexcel D3 = "coef"
putexcel E3 = "%"
putexcel F2 = "$\rho = 0.8$"
putexcel F3 = "coef"
putexcel G3 = "%"

putexcel H1 = "T = 100"
putexcel H2 = "$\rho = 0.2$"
putexcel H3 = "coef"
putexcel I3 = "%"
putexcel J2 = "$\rho = 0.5$"
putexcel J3 = "coef"
putexcel K3 = "%"
putexcel L2 = "$\rho = 0.8$"
putexcel L3 = "coef"
putexcel M3 = "%"

putexcel N1 = "T = 200"
putexcel N2 = "$\rho = 0.2$"
putexcel N3 = "coef"
putexcel O3 = "%"
putexcel P2 = "$\rho = 0.5$"
putexcel P3 = "coef"
putexcel Q3 = "%"
putexcel R2 = "$\rho = 0.8$"
putexcel R3 = "coef"
putexcel S3 = "%"



*************************
**T=50; x1, rho=.2, .5, .8
*************************
set seed 457090451

******
*rho=.2
*****
*program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 50
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.2*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

*Simulate the program "combined" N times and save the betas and standard errors.
*Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
*can correctly identify TRUE relationships
simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************


*Generate t-statistic for each simulated regression and export results into excel

*export results from ADL rho = 0.2 to excel file
**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel B4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel C4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel B5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel C5 = `perc'

**$\hat{\beta}_2$
*ceof
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel B6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel C6 = `perc'


******
*\rho=.5
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 50
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.5*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************

*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.5 to excel file
**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel D4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel E4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel D5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel E5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel D6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel E6 = `perc'

******
*rho=.8
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 50
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.8*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************

*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.8 to excel file
**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel F4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel G4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel F5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel G5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel F6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel G6 = `perc'


*************************
**T=100; x1, rho=.2, .5, .8
*************************

******
*rho=.2
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 100
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.2*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

*Simulate the program "combined" N times and save the betas and standard errors.
*Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
*can correctly identify TRUE relationships
simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq


*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************

*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.2 to excel file
**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel H4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel I4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel H5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel I5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel H6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel I6 = `perc'


******
*rho=.5
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 100
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.5*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************

*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.5 to excel file
**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel J4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel K4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel J5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel K5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel J6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel K6 = `perc'


******
*rho=.8
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 100
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.8*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************

*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.8 to excel file

**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel L4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel M4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel L5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel M5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel L6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel M6 = `perc'



*************************
**T=200; x1, rho=.2, .5, .8
*************************

******
*rho=.2
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 200
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.2*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

*Simulate the program "combined" N times and save the betas and standard errors.
*Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
*can correctly identify TRUE relationships
simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************
*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.2 to excel file

**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel N4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel O4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel N5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel O5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel N6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel O6 = `perc'



******
*rho=.5
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 200
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.5*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
**%
*****************
*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.2 to excel file

**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel P4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/1000
di `perc' 
putexcel Q4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel P5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/1000
di `perc'
putexcel Q5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel P6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/1000
di `perc'
putexcel Q6 = `perc'



******
*rho=.8
*****
program drop combined_noq
program define combined_noq, rclass
drop _all
set obs 200
gen t = _n
*gen stationary time series (x1)
gen e1=invnorm(uniform())
gen x1=e1 if t==1
replace x1=.8*x1[_n-1] + e1 if t>1
*gen integrated time series (x2)
gen u=invnorm(uniform())
gen x2=u if t==1
replace x2=x2[_n-1] + u if t>1
*gen combined times series (z) that his a function of x1 and x2
gen q=invnorm(uniform())
gen z = x1 + x2
tsset t
reg z l.z x1 l.x1
	estat bgodfrey
	mat P = r(p)
	return scalar pvalue_bg = P[1,1] 
end

simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq

*Percent of simulations which a Breusch-Godfrey test rejects teh null of 
*no serial correlation
sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

*****************
**coef
*****************
sum _sim_1 _b_x1 _sim_3

*****************
** %
*****************
*Generate t-statistic for each simulated regression and export results into excel
*export results from ADL rho = 0.2 to excel file

**$\hat{\alpha}_1$
*coef
quietly sum _sim_1
di %6.4f `r(mean)'
putexcel R4 = `r(mean)', nformat("#0.####")
*%
gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
quietly sum tstat_ly
local perc = 100*`r(N)'/2000
di `perc' 
putexcel S4 = `perc'

**$\hat{\beta}_1$
*coef
quietly sum _b_x1
di %6.4f `r(mean)'
putexcel R5 = `r(mean)', nformat("#0.####")
*%
gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
quietly sum tstat_x1
local perc = 100*`r(N)'/2000
di `perc'
putexcel S5 = `perc'

**$\hat{\beta}_2$
*coef
quietly sum _sim_3
di %6.4f `r(mean)'
putexcel R6 = `r(mean)', nformat("#0.####")
*%
gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
quietly sum tstat_lx1
local perc = 100*`r(N)'/2000
di `perc'
putexcel S6 = `perc'

